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We experimentally study the effects of coupling one-dimensional Many-Body Localized (MBL) 
systems with identical disorder. Using a gas of ultracold fermions in an optical lattice, we artifically 
prepare an initial charge density wave in an array of ID tubes with quasi-random onsite disorder and 
monitor the subsequent dynamics over several thousand tunneling times. We find a strikingly dif¬ 
ferent behavior between MBL and Anderson Localization. While the non-interacting Anderson case 
remains localized, in the interacting case any coupling between the tubes leads to a delocalization 
of the entire system. 


Introduction. —Many-Body Localization (MBL) marks 
a new paradigm in condensed matter and statistical 
physics. It describes an insulating phase in which a dis¬ 
ordered, interacting many-body quantum system fails to 
act as its own heat bath [1-5]. In isolation, these sys¬ 
tems will never achieve local thermal equilibrium and 
conventional statistical physics approaches break down. 
Unlike other insulating phases, MBL is not limited to 
ground states but can even occur in all exited states 
of a disordered many-body system [6-10]. A dynami¬ 
cal phase transition separates the MBL phase from con¬ 
ventional ergodic phases [11, 12], in which the isolated 
system thermalizes. In these ergodic phases, any initial 
quantum information becomes rapidly diluted in the ex¬ 
ponentially large Hilbert space, leading to decoherence. 
In contrast, in the localized insulating phase quantum 
information can persist locally for an infinite amount of 
time [5]. This could potentially render quantum infor¬ 
mation devices less susceptible to noise and disorder. For 
many decades it remained unclear whether such a local¬ 
ized phase could persist in a many-body system beyond 
the non-interacting limit of an Anderson insulator [13]. 
Today, both theory and experiment have shown evidence 
for the existence of an MBL phase in interacting ID sys¬ 
tems [5-7, 14-16]. Nonetheless, many fundamental ques¬ 
tions regarding this phase and the associated phase tran¬ 
sition as well as its extension to higher dimensions [17] 
remain open, making it a highly active topic of current 
theoretical and experimental research. 

One crucial requirement for the existence of a MBL 
phase is that no coupling to any external heat bath or 
bath-like structure exists. Any such coupling will even¬ 
tually thermalize the system and ultimately destroy the 
non-thermalizing MBL phase [5, 18-22], Since any exper¬ 
imental system will inevitably be coupled—albeit poten¬ 
tially very weakly—to an environment, it is of critical im¬ 
portance to quantitatively understand the effect of such 
a coupling. Furthermore, studying the effects of weak 
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FIG. 1. Coupling identical MBL systems: A charge density 
wave (CDW) with atoms only occupying even sites (e) is pre¬ 
pared in each of the identically disordered ID tubes along the 
longitudinal ( x ) direction, with hopping J, on-site interaction 
energy U and disorder strength A. Red and blue spheres in¬ 
dicate a typical distribution of |t) and |4) atoms. We monitor 
the time evolution of such a state for different inter-tube cou¬ 
pling strengths J±, that is different hopping amplitudes along 
the transverse (y) direction. 


couplings can help to experimentally identify an MBL 
phase and distinguish it from glasses or non-interacting 
Anderson localized phases. 

Ultracold quantum gases in optical lattices form an 
ideal system to investigate these questions, as they are 
almost perfectly isolated from the external world and 
are highly controllable. Earlier experiments have in¬ 
vestigated the interplay between disorder and interac¬ 
tions in the ground state of an isolated systems of ul¬ 
tracold bosons in optical lattices [23-25] and studied the 
influence of disorder on transport properties for lattice 
fermions [26]. Recently, we were able to show that for 
a wide range of energy-densities and interactions, an 
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MBL phase exists in ID Hubbard type chains with quasi¬ 
random disorder [14]. 

In this work, we experimentally study the effects of 
coupling identically disordered ID MBL systems to each 
other. In particular, the disordered Hamiltonian is identi¬ 
cal for all ID tubes, but the initial configuration of atoms 
differs between the tubes (Fig. 1). We find that the cou¬ 
pled systems can collectively serve as a bath for each 
other, i.e. coupling localized systems can result in delo¬ 
calizing all of them. 

Experiment. —Our experiments start with a two com¬ 
ponent Fermi gas of 40 K atoms in an equal mixture of 
the two lowest hyperhne states |F, m F ) = |§,— §) = |i) 
and ||,— |) = |t) with a total atom number of about 
110-150 x 10 3 atoms. In the initial dipole trap, the 
atoms are at a temperature of 0.19(2) Tp, where T F is 
the Fermi temperature. We load the Fermi gas into the 
lowest band of a deep, three-dimensional simple cubic op¬ 
tical lattice, where tunneling can be neglected. Along the 
longitudinal (x) direction, we then add a second (short) 
lattice (wavelength A s =532nm) to the initial (long) lat¬ 
tice (A;=1064nm). By controlling the phase of the short 
lattice during loading, we prepare a period-two ’Charge- 
Density-Wave’ (CDW), where only even sites are occu¬ 
pied in the 20 E r deep short lattice. Here, the recoil 
energy is denoted by E r = h 2 /2mA 2 , where h is Planck’s 
constant, A is the respective lattice wavelength and m is 
the atomic mass. The orthogonal lattices along y and z 
with a wavelength of A_l= 738 nm are initially ramped up 
to 45 E r , creating an array of (almost) isolated ID tubes. 
During lattice loading the interactions are kept strongly 
repulsive at a scattering length of a = 142 ao, where ao is 
the Bohr radius, by employing a Feshbach resonance cen¬ 
tered at 202.1 G [27]. This results in a doublon fraction, 
that is the fraction of atoms on doubly occupied lattice 
sites, of < 10%. 

After the preparation of the CDW in the deep lattices, 
the desired interactions for the ensuing evolution are set. 
Additionally, an incommensurate lattice of wavelength 
Ad=738 nm is superimposed along the x direction to cre¬ 
ate quasi-random onsite disorder along the longitudinal 
direction. The system size is approximately 200 sites in 
the longitudinal and 120 sites in the transverse direc¬ 
tion, and the central longitudinal tubes contains about 
90 atoms [28]. After this preparation, the long lattice is 
quickly ramped to zero, the transverse y lattice is ramped 
to its final value, which controls the transverse coupling 
Ji_, and the short lattice is reduced to 8 E r . This last 
ramp enables tunneling along the tube and thereby initi¬ 
ates the dynamics. After a variable evolution time, we ex¬ 
tract the imbalance between atoms on even and odd sites 
I = (N e — N 0 )/(N e + N a ) [14]. Here N e and N 0 denote 
the population of even and odd sites respectively along 
the longitudinal direction and are extracted by mapping 
them to different bands of the superlattice [29]. 

The imbalance provides a measure of ergodicity break- 



FIG. 2. Time evolution of a charge-density wave (CDW): An 
initially prepared ID CDW evolves in a disordered system 
with disorder strength A = 5 J along one axis. We measure 
the imbalance I after a given evolution time in 1D + (circles), 
where J± < 1CC 3 J and in 2D (diamonds), where J± = J. 
Each data point is the average of 6 disorder phase realiza¬ 
tions, with error bars showing the standard error of the mean. 
Solid lines are fits [28] from which we extract the imbalance 
lifetimes. Shown in gray is an exact diagonalization (ED) 
calculation for J± = 0, A = 5 J and U = 0 [28]. 

ing: It quickly decays to zero under any ergodic dynamics 
and any non-zero imbalance persisting at long times sig¬ 
nifies a memory of the initial state and directly indicates 
localization [14, 22, 30]. 

Model. —Our system consists of an array of ID tubes 
with identical quasi-random on-site disorder along the 
longitudinal direction. Each tube can be described by 
the Aubry-Andre Model [31] with interactions [32], as de¬ 
picted in Fig. 1. A finite hopping amplitude J± along the 
transverse (y) direction introduces a coupling between 
adjacent tubes. The Hamiltonian of the system is given 
by 

H = — J .o-Gj.j + h.c.) 

~ j± 1,oAj> + h' C ') 

i,j, cr 

cos(27 t/3i + (j>)hi j >cr + U ^ nqytfijjq, 

i,j,cr i,j 

(1) 

where J ~ h x 500 Hz is the tunneling matrix element 
between neighboring sites along a tube and J± denotes 
the transverse ( y ) hopping between the tubes. The cre¬ 
ation (annihilation) operator for a fermion in spin state 
a £ {t)l} on site i in tube j is c\ j a (ci,j,a) and the lo¬ 
cal number operator is given by ni,j >a = c| • a Ci,j,a- The 
quasi-random on-site disorder is characterized by the dis- 
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FIG. 3. Imbalance lifetimes versus interactions (U): Imbal¬ 
ance lifetimes at A = 5 J for 1D + and 2D. We note that the 
1D + case differs crucially from the ideal isolated ID case. The 
lifetimes were extracted from fits to time traces such as in Fig. 
2. Error bars denote fit uncertainty [28]. The gray shaded 
area indicates the range of measured atom number lifetimes, 
while the lines are guides to the eye. 


order amplitude A, the incommensurable wavelength ra¬ 
tio /3 = A s /A d and the relative phase (j). Finally, the 
on-site interaction energy is given by U . 

In the limit Jj_ —> 0, the system decouples into many 
one-dimensional tubes, which show many-body localiza¬ 
tion [14]. For our experiment, the accessible limits of 
almost zero inter-tube coupling J± < 10~ 3 J and equal 
coupling J± = J are termed the 1D + and 2D cases, re¬ 
spectively. Note that the experimentally achievable limit 
of the 1D + case is not the same as the ideal theoretical 
ID case. A small but non-zero coupling remains for any 
finite transverse lattice depth and can always affect the 
dynamics at very long time scales. 

Results. —We monitor the time evolution of the imbal¬ 
ance at disorder amplitude A = 5 J for various interac¬ 
tion strengths U in both the 1D + and 2D cases. This 
disorder strength is deep in the MBL regime for isolated 
ID tubes. Fig. 2 shows exemplary time traces at U = 0 
and U = —6 J, with all times given in units of the longi¬ 
tudinal tunneling time r = h/ J . We start with an out- 
of-equilibrium density wave with an initial imbalance of 
I(t = 0) = 0.91 ± 0.03 and observe a fast initial decrease 
of the imbalance up to approximately one tunneling time. 
This decrease is similar in all cases and corresponds to 
an initial relaxation in the longitudinal direction. 

As shown in Fig. 2, in both non-interacting cases, the 
initial decrease is followed by highly damped oscillations 
around a plateau at finite imbalance, closely matching 
the expected steady state for this Anderson localized sys¬ 
tem [14, 28]. At very long times (» 100 r), the curves 
start to deviate from this plateau and exhibit a slow de¬ 
cay. The corresponding lifetime is extracted by fitting 


FIG. 4. Imbalance lifetimes in the 1D-2D crossover: The 
coupling strength J± between adjacent tubes is varied con¬ 
tinuously for A = 5 J at four different interactions U. Here 
J±/J < 10 -3 and J±/J = 10° correspond to the 1D + and 
the 2D cases respectively. The lack of saturation as J± —> 0 
indicates that the residual inter-tube coupling still limits the 
imbalance lifetime in the 1D + case. Solid lines denote power- 
law fits oc J[5[ for small Jj_, with fitted exponents k(U/J) of 
k( 0) = 0.00(4), k(- 2) = -0.09(2), fc(-6) = -0.30(1) and 
k(— 20) = —0.16(3). We note that in principle the tunneling 
along the z direction becomes sizable for the smallest J±, as 
Ji/J~ 1(T 3 . 

the imbalance traces to a damped sinusoid, which models 
the initial fast relaxation, multiplied by a stretched ex¬ 
ponential to capture the slow decay [28]. This stretched 
exponential is of the form e, where t is the evolution 
time, r is the decay rate and /3 is the stretching exponent. 
This model fits all time traces consistently better than a 
simple exponential decay and has also been observed in 
glasses, disordered materials and polymers [33-35]. We 
define our imbalance lifetime Tl/e = 1/r. 

In the absence of interactions, the observed long time 
dynamics is dominated by classical noise, photon scatter¬ 
ing from lattice beams and other technical imperfections. 
These processes couple the system to the environment 
and over time delocalize it [36, 37]. In addition, these 
experimental imperfections also give rise to an atom num¬ 
ber decay, which limits the lifetime of atoms in the lattice 
to 0.5 —1.1 x 10 4 r [28]. In the absence of interactions, the 
observed dynamics in 1D + and 2D are almost the same, 
since the disorder potential is identical in all tubes [28]. 
Therefore the 2D system is separable, and the longitudi¬ 
nal and transverse directions are decoupled. Note that 
we do not expect this to hold if the disorder were differ¬ 
ent in different tubes [38]. Adding interactions, however, 
breaks the separability of the system and the transverse 
dynamics can now affect the imbalance along the longi¬ 
tudinal direction. Since there is no disorder along the 
transverse direction, particles are free to move along this 
direction. In the interacting case, this couples the orig- 
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inally localized tubes such that they collectively act as 
a bath for each other and thereby delocalize the entire 
system. As a consequence, the interacting 2D trace in 
Fig. 2 displays no plateau but instead shows a fast de¬ 
cay. To understand this behavior further, we map out the 
imbalance lifetimes in the 1D + and 2D cases for various 
interaction strengths (Fig. 3). We find that in 2D, even 
small interactions are sufficient to dramatically reduce 
the imbalance lifetime to less than 100 r. 

Although MBL is expected to be stable in the isolated 
ID case, the 1D + data in Fig. 3 shows qualitatively sim¬ 
ilar behavior to the 2D case, but with a much weaker de¬ 
crease of imbalance lifetimes with increasing interactions. 
This suggests that the small but non-zero inter-tube cou¬ 
pling also limits the lifetime in the 1D + case. The gray 
shaded region in Fig. 3 marks the range of observed atom 
number lifetimes, which approximately coincides with the 
non-interacting imbalance lifetimes, indicating the rele¬ 
vance of technical imperfections on this timescale. 

1D-2D Crossover.- In order to directly test the effect 
of residual couplings, we vary the strength of the inter¬ 
tube coupling between the 2D and the 1D + limits for 
four different interactions (Fig. 4). We observe increas¬ 
ing lifetimes for decreasing coupling strengths J_l in all 
interacting cases. For small but finite coupling strengths, 
we observe a linear trend on a log-log scale (Fig. 4), sug¬ 
gesting a power law dependence. For strong inter-tube 
coupling (J±/J > 0.1), there is a crossover to a faster de¬ 
cay. The fitted exponents are surprisingly small (\k\ < 1) 
and depend non-monotonically on the interactions. In 
the non-interacting limit C/ —>• 0, the lifetimes are on the 
order of the atom number lifetime and become indepen¬ 
dent of the transverse coupling due to the separability 
of the problem. This highlights the striking difference 
between MBL and Anderson localization. 

Extrapolating towards the limit of the true ID case 
(Jj_ —> 0), we would expect the lifetimes of the interact¬ 
ing system to saturate once the inter-tube coupling is no 
longer the dominant decay mechanism. However, in the 
experimentally accessible regime we cannot observe any 
signs of saturation, strongly suggesting that the inter¬ 
acting lifetimes are still limited by the non-zero residual 
transverse coupling to the neighboring tubes, even in the 
1D + case. 

Constant evolution time. —The interaction dependence 
is also visible in the imbalance value measured after a 
fixed evolution time of around 40 t, as shown in Fig. 5. 
We observe a substantial decrease of the 2D imbalance 
with small interactions, which is consistent with the 
sharply decreasing lifetimes. Additionally, Fig. 5 shows 
that the dynamics are symmetric around U = 0, which 
is expected due to a dynamical symmetry of the Fermi- 
Hubbard Model [39]. Interestingly, we observe an in¬ 
creasing imbalance in the 2D data for strong interactions. 
We checked that this increase is not due to doublons 
by removing any residual doublons with a pulse of near- 



FIG. 5. Imbalance vs. interactions ( U ) at fixed times: Data 
taken at A = 5 J averaged over three different times (38 r — 
41 r) and 4 disorder phases 4>. In the 1D + case (circles), the 
imbalance corresponds to the stationary value [14], whereas in 
the 2D case (diamonds) it is indicative of the decay lifetime. 
Solid lines are guides to the eye. The 2D case is limited to 
U/J < 10 due to details of the used Feshbach resonance [27]. 

resonant light [40] prior to the evolution. This increase 
in lifetime might be due to the reduced phase space for 
scattering in the hard-core limit (U J). A recent 

theoretical study using a cluster expansion method on a 
smaller system (8x8 sites) observed a similar trend [41]. 

At 40 r, the 1D + case is dominated by the ID plateau 
value and shows the characteristic ‘W’ shape of the 
steady state imbalance of the MBL system [14] , with lit¬ 
tle influence from the (much longer) lifetimes. As shown 
in the supplementary material [28], we additionally find 
that the imbalance lifetime in all cases increases strongly 
for larger disorder. 

Conclusion .—We have studied the stability of many- 
body localized ID systems under an inter-tube coupling. 
We found that even weak couplings have a delocalizing ef¬ 
fect on the MBL phase, while leaving the non-interacting 
Anderson limit unchanged. This highlights the differ¬ 
ences between these two regimes and shows the princi¬ 
pal fragility of MBL with respect to coupling to any ex¬ 
ternal heat bath [5, 18-21], Furthermore, we have not 
observed any saturation in imbalance lifetimes even for 
the smallest couplings we could attain, indicating that 
this inter-tube coupling is the dominant decay mecha¬ 
nism in this MBL experiment. Nonetheless, for strong 
disorders the achieved lifetimes already exceed the life¬ 
times of typical many-body states, such as e.g. superfluid 
states in optical lattices. This demonstrates the stabil¬ 
ity of MBL with respect to other experimental imperfec¬ 
tions and increases the prospects of realizing localization 
protected order [42-44] and applications in quantum- 
information [5]. 

An important next step will be to extend this study to 
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the ‘true’ 2D case with disorder along both directions. In 
addition, future experiments should also be able to ad¬ 
dress the question of the stability of MBL under external 
influences, such as photon scattering [36, 45] and time- 
dependent modulations [46, 47]. Finally, it would also be 
interesting to search for MBL in bosonic systems [48] and 
using other observables such as the growth of entangle¬ 
ment entropy and response to generalized interferometric 
probes [49] when a bath is added to an otherwise perfectly 
MBL system [22]. 
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and A. Pal. We acknowledge financial support by the Eu¬ 
ropean Commision (UQUAM, AQuS) and the Nanosys¬ 
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Supplementary Material 

Imbalance Lifetimes vs Disorder Strength. —Fig. SI 
shows the dependence of imbalance lifetimes on the disor¬ 
der strength A in the 1D + and the 2D case for U = —6 J. 
In both cases, we observe exponentially increasing life¬ 
times with larger A. We attribute this to a reduced effec¬ 
tiveness of any bath-like perturbation for large disorder 
strengths. At long lifetimes this behavior saturates, as 
the imbalance lifetimes approach the same order of mag¬ 
nitude as the atom loss lifetimes. Since all time traces in 
the paper are measured up to 3000 r, lifetimes exceeding 
5000 r may be affected by systematic errors due to the 
fitting procedure becoming unreliable. 

To capture the behavior with A we fit the data with 
the following empirical function: 

T v (S01) 

where Tf, represents the background imbalance loss 
due to noise, while a and b are free parameters. These 
fits are represented as solid lines in Fig. SI. The fit 
routine is performed in logarithmic units. 

Non-Interacting Lifetime Difference between 1D + and 
2D. The difference between the non-interacting 1D + 
and the 2D data is most likely due to technical imperfec¬ 
tions. Firstly, the interactions U may not be exactly zero 
for all times as it relies on the stability and the absolute 
calibration of the magnetic field. Secondly, the lattice 
beams have gaussian envelopes. This implies that both 
longitudinal hopping and disorder depend slightly on the 
transverse position. Thirdly, the relative phase of the 
disorder might be slightly different between neighboring 
tubes due to slight misalignment between the disorder 
beam and the main lattice. These imperfections break 
the separability of the problem. 

Atom Number Lifetimes. - In order to determine the 
atom number lifetime, we fit an exponential decay to the 
measured atom number time trace for every imbalance 
dataset. Five typical curves are shown in Fig. S2. The 
lifetimes are similar across the full investigated parameter 
range and we see no strong trend with either A or U. 
Extracted lifetimes from the exponential fits lie in the 
range of 0.5 — 1.1 x 10 4 r. 

As the experimental lifetime of atoms in the pure 
dipole trap is over 10 5 r, the detected atom loss in the 
lattice must be dominated by off-resonant photon scatter¬ 
ing of lattice photons as well as frequency and amplitude 
noise in the optical fields. This is consistent with the 
extracted atom number lifetimes being slightly higher in 
the 2D cases as compared to the 1D + cases, as the to¬ 
tal optical power in the lattice beams is lower in the 2D 
cases. The loss channels formed by such processes will 
most likely involve excitations of atoms to higher bands. 



FIG. SI. Imbalance lifetimes as a function of disorder 
strength: The lifetimes were extracted by fits to individual 
time traces taken up to 3000 r, as in the main text. Er¬ 
ror bars indicate uncertainty in the fitting procedure. The 
gray shaded area indicates the range of typical atom num¬ 
ber lifetimes. The lifetimes for the non-interacting cases at 
larger disorder strengths exceed the atom number lifetimes 
and hence cannot be reliably extracted. Solid lines show the 
results of saturated exponential fits as described in the text. 

We note that, due to the long observed atom number 
lifetimes, these loss channels are unlikely to be relevant 
for the fast imbalance decay at the observed timescales in 
2D: Previous DMRG simulations in ID [14] agreed well 
with the experimental measurements in the 1D + case un¬ 
til « 100 r. Hence, we do not expect higher bands to 
play a significant role up to at least this time. Further¬ 
more, the atom number lifetimes in the 2D cases are even 
slightly longer than in the 1D + cases. Besides, at a rela¬ 
tively short time of 40 r, as presented in Fig. 5, we already 
find a strong effect of the inter-tube coupling, despite the 
fact that at these time scales there is no discernible atom 
loss. 

Extracting Imbalance Lifetimes. — The observed time 
evolutions of the imbalance X displays some common 
characteristic features. The imbalance starts close to 
unity, followed by a couple of oscillations with a period 
on the order of r, which quickly settle to a non-zero 
plateau value, before decaying to zero in the long-term 
limit. To quantify this long term decay we fit the sum 
of a damped sinusoidal and an offset multiplied with an 
overall stretched exponential decay. This fitting function 
is given by: 

X = {Ae~ t ^ tl cos{vt) + / st )e _(r ^ /3 (S02) 

where A is the amplitude of the sinusoid with oscillation 
frequency v ~ 1/r, oscillation damping time t\ ~ r, 
plateau value I s t, long-term decay rate T and stretching 
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FIG. S2. Normalized atom loss time traces: Characteristic 
traces showing the loss of atoms for different couplings J_l, 
disorders A and interactions U. Every data point is averaged 
over 6 disorder realizations <j>. Each trace was individually 
normalized to its maximum value after averaging. Error bars 
show the error of the mean. Lines are fit with exponentials of 
the form Ae _t / T “. 


coefficient ft. From this fit function we obtain the 
lifetimes TV e = 1/r, as shown in Figs. 3, 4 and SI. This 
lifetime would diverge in the case of a perfectly localized 
system without external noise and strictly vanishing 
bath couplings. The error bars in T\i e throughout 
the paper are fit uncertainties. We fit a stretched 
exponential instead of a pure exponential as it captures 
the measured time traces better. This is illustrated in 
Fig. S3. Note that the choice of the stretched or the 
pure exponential does not qualitatively change any the 
observed behavior. 

Exact Diagonalization .—The ED calculations are 
adapted from our previous work. We simulate the 
Anderson localization limit by performing Exact Diago¬ 
nalization (ED) on a single, non-interacting particle and 
account for our experimental setup by averaging over 
fifty phases (equally spaced between [0, 2tv\), the dipole 
trap, the inhomogeneity of the lattice and the atomic 
cloud shape. Full details of the ED can be found in 
Ref. [14]. 

General Sequence .—The experiment produces an ul¬ 
tracold gas of fermionic Potassium-40 ( 40 K) atoms by 
sympathetically cooling 40 K with bosonic Rubidium-87 
( 87 Rb) in a plugged quadrupole trap followed by an opti¬ 
cal dipole trap. Reducing the dipole trap depth lower 
than a threshold value completely removes 87 Rb due 
to its higher mass, such that only 40 K remains in the 
trap. We further evaporate 40 K in an equal mixture 
of the two lowest hyperfine states of the A S 1 / 2 manifold 
(I F,m F ) = |9/2, —9/2) = |4) and |9/2, -7/2) = |f» to a 


final temperature of T/Tp = 0.19(0.02), where Tp is the 
Fermi temperature. Interactions between the two states 
can be tuned via an s-wave Feshbacli resonance centered 
at 202.1G [27]. The scattering length a is set to 140 a 0 
during lattice loading to avoid double occupancies. Here, 
ao denotes the Bohr radius. The number of double occu¬ 
pancies are characterized by converting them into Fesh- 
bach molecules via crossing the Feshbach resonance and 
subsequently imaging the remaining single atoms [50]. 

After loading into the deep lattices (as described in the 
main text), we then set the scattering length to control 
the desired interactions U in the lattice for the following 
evolution time. At the end of this preparation stage, 
we ramp down the x-short lattice from 20 E r to 8 E r in 
10 fj,s and ramp down the rr-long lattice to 0 E r in the same 
time. The disorder lattice is simultaneously ramped up to 
the desired disorder strength. The orthogonal y-lattice is 
ramped to the final value with an additional time of 90 fjs. 
All of these lattice ramp times are short compared to a 
tunneling time r. The system is then allowed to evolve for 
various evolution times. For detection, the short lattice, 
the long lattice and the orthogonal lattices are ramped 
high again to inhibit hopping and freeze all occupations. 
Finally, we employ the bandmapping technique to obtain 
the number of atoms on the even and odd sites [29]. All 
bandmapping images are taken with 8 ms time of flight 
and are imaged along the y-axis (i.e. orthogonal to the 
superlattice axis). 

The cloud size in the lattice is measured via in-situ 
imaging. It can be described by a gaussian, with 
standard deviation a of ~ 14 //m and ~ 11 /tm in the x 
and y directions respectively. To estimate the number 
of sites in the central lattice tube, we take a value which 
corresponds to 4cr, i.e. the range from — 2er to +2 a. 

Lattice Parameters. - All lattice potentials result from 
retroreflected laser beams with gaussian intensity profiles 
and 1/e 2 radii (waists) of 150/xm that are centered at 
the position of the atoms. The lattice beams along the 
y and z direction and the disorder lattice light originate 
from the same Coherent MBR-110 Ti:Sa laser operated 
at 738 nm, which is locked to its internal reference cavity. 
In order to eliminate cross-interference between beams 
along different axes, separate beams have orthogonal 
polarizations and are shifted to different frequencies 
using acousto-optic modulators. The relative phase <j> of 
the disorder potential is controlled by slightly changing 
the frequency of this laser. The 1064 nm laser used for 
the long wavelength superlattice along the x-direction is 
locked to a dedicated reference cavity that is acoustically 
isolated and thermally stabilized to 65 kHz linewidth 
over 100 ms. A 532 nm Coherant Verdi laser is used as 
the short lattice and is locked to frequency doubled light 
from the 1064 nm laser via an offset lock. 

Disorder Strength. —The disorder strength A depends 
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FIG. S3. We reproduce the time traces from Fig. 2 in the 
main paper with both the stretched exponential (solid lines, 
eq. S02) and a pure exponential fit (dotted lines) demonstrat¬ 
ing that the former describes the measured data better. We 
find this behavior to hold for all the time traces we have mea¬ 
sured. 



FIG. S4. Dynamical Symmetry of the Hubbard Model: Time 
traces showing the evolution of the imbalance X in 2D for 
A = 5 J and U = ±10 J. The data is identical within ex¬ 
perimental scatter, showing that the decay of the imbalance 
does not depend on the sign of U . Each point is averaged 
over six disorder realizations q i. Error bars denote the error 
of the mean. Solid lines are fitted stretched exponentials as 
per eq. S02. 

on the lattice depth of the 532nm laser (main lattice), 
the lattice depth of the disorder lattice and the ratio of 
their wavelengths f3 ~ 532/738. We use 8 -E r ,532 as the 
main lattice depth, where -EV -,532 = h x 17.64 kHz is the 
recoil energy in the 532 nm lattice. The disorder strength 
is then given by A /J = 6.67 s<j, where Sd, is the depth 
of the disorder lattice in units of the 738 nm recoil energy. 

Dynamical U <± —U Symmetry .—We load the atoms 
in the lattice at magnetic fields above the Feshbach 



FIG. S5. Fitted Power Law Exponents for 1D-2D Crossover: 
Shown are the power law exponents of the 1D-2D crossover 
as a function of interaction strength. Error bars denote the 
fit uncertainty. 

resonance at 202.1 G [27]. On this side of the resonance, 
the maximum positive scattering length that can be 
obtained with our experiment is a ss 142 ao- This 
scattering length corresponds to an interaction strength 
of U sa 12 J in the 2D case. However, there is no such 
restriction in the case of negative interactions. Hence we 
measure most of the data with attractive interactions 
(i.e. U negative). From a dynamical symmetry of the 
Hubbard model [39], we expect the dynamics to be 
identical for ±U, as was observed to be the case for 
the short time behavior in the 1D + case [14]. However, 
to investigate whether the symmetry also holds for the 
long term dynamics, we measure the time traces in 
the 2D case for U = ±10 J. This data is shown in 
Fig. S4 and indicates that the relaxation dynamics are 
indeed symmetric about U = 0. Since fluctuations of 
the magnetic field would cause much smaller changes 
of the scattering length in the U = ±10 J case than in 
the U = —10 J case, the observed U <± —U symmetry 
indicates that neither field fluctuations nor the proximity 
to the Feshbach resonance contribute significantly to the 
decay of imbalance in the interacting 2D cases. 

Doublon Dependence. —While all the data presented in 
this paper is taken with a rather small doublon fraction 
of ~ 8 %, we also briefly investigate the effect of doublons 
on the imbalance lifetimes. We measure the lifetimes 
in the 1D + and the 2D case at A = 5 J, U = —6J 
for doublon fractions of 8 % and 50 % as well as for 
U = — 20 J with 0% and 8% doublons. We can vary 
the doublon fraction between ss 8% to « 50% by 
varying the scattering length during lattice loading 
between ai oa( i = —140 ao to ±140ao- In addition, we can 
reduce the doublon fraction further to almost zero by 
applying a near resonant blast pulse in the deep lattice. 
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singly occupied sites unaffected. For all of the above 
cases, we observe less than 5 % change in the imbalance 
lifetime, with slightly longer imbalance lifetimes for 
higher doublon fractions. 

Power-Law Fits for 1D-2D Crossover. —The solid 
lines shown in Fig. 4 are fits to Ti/ e = b(J±/J) k , where 
amplitude b and the power law exponents k are free 
fit parameters. The fitted exponents are plotted in 
Fig. S5 as a function of interactions. Note that the 
exponents for all the interacting cases are negative with 
absolute values much smaller than unity. The exponents 
depend non-monotonously on the interaction strength, 
resembling the shape of the 2D curve in Fig. 5. 


FIG. S6. Variation of the stretching exponents /3 in the 1D- 
2D Crossover: Shown are the fitted stretching exponents (5 
in the 1D-2D crossover as a function of interaction strength. 
Error bars denote the fit uncertainty. 


This removes the leftover doublons while leaving the 


Variation of the stretching exponents in the 1D-2D 
crossover. — The fitted stretched exponents /3 are shown 
in Fig. S6 and typically lie in the range [0.5,1]. We do 
not observe a systematic dependence of the stretching 
exponent on the transverse coupling J±. However, the 
exponents seem to be systematically larger in the case 
of U = -20 J. 












